Conformational Heterogeneity and Interchain Percolation Revealed in an Amorphous Conjugated Polymer

Conjugated polymers are employed in a variety of application areas due to their bright fluorescence and strong biocompatibility. However, understanding the structure of amorphous conjugated polymers on the nanoscale is extremely challenging compared to their related crystalline phases. Using a bespoke classical force field, we study amorphous poly(9,9-di-n-octylfluorene-alt-benzothiadiazole) (F8BT) with molecular dynamics simulations to investigate the role that its nanoscale structure plays in controlling its emergent (and all-important) optical properties. Notably, we show that a giant percolating cluster exists within amorphous F8BT, which has ramifications in understanding the nature of interchain species that drive the quantum yield reduction and bathochromic shift observed in conjugated polymer-based devices and nanostructures. We also show that distinct conformations can be unravelled from within the disordered structure of amorphous F8BT using a two-stage machine learning protocol, highlighting a link between molecular conformation and ring stacking propensity. This work provides predictive understanding by which to enhance the optical properties of next-generation conjugated polymer-based devices and materials by rational, simulation-led design principles.

Solvation Mechanisms of F8BT in Water and THF Figure S3 shows radial distribution functions (g(r)) and associated coordination numbers that describe the solvation mechanisms of F8BT in water ( Figure S3(a)) and THF ( Figure   S3(b)). F8BT is insoluble in aqueous solution. In water, the benzothiadiazole ring is more strongly hydrated than the fluorene moiety. The highest level of solvation in the benzothiadiazole ring is seen for the S atom. Similarly, we observe that the C D and C E atoms are the most hydrated fluorene ring atoms, as these atoms are those least shielded by the octyl side chains. We observe preferential solvation of specific benzothiadiazole and fluorene atoms in water, while for F8BT in THF ( Figure S3  Similar relative behaviour is also observed in the solvation mechanisms of the octyl side chains in water and THF ( Figure S4). The coordination number of THF around the alkyl chain increases monotonically moving down the chain from the fluorene ring (with the exception of the terminal methyl group, which is expected to be due to the slight difference in steric hindrance conferred by the methyl H atoms). Conversely, for the case of F8BT in water, we note that preferential hydration is experienced by C3 compared to atoms we might expect to be found further into aqueous solution (C4, C5, and C6). This indicates that hydrophobic collapse of the side chains reduces solvent interactions with the side chains.
(a)  Figure S4: F8BT side chain solvation in (a) water and (b) THF. A cutoff value of 7.5Å was used to calculate the coordination numbers shown as bar charts. RDFs are vertically shifted by 1 for clarity.

Quantifying Polymer Deformation and Relaxation
The following quantity provides a simple metric by which to compare the dynamics of polymer chains in different conditions (melt phase and in solution) and at different temperatures. 1 The root-mean squared deviation (RMSD) of a polymer chain given a lag time, τ , is defined Now we define Γ(τ ) ≡ ⟨RMSD (r(t), r(t + τ ))⟩, the RMSD as a function of the lag time, averaged over all molecules and different start times t. Assuming that Γ → Γ max as τ → ∞, we model the relaxation process using a simple first order ordinary differential equation Integrating the ODE and applying the initial conditions Γ = 0 when τ = 0 yields We evaluate A = Γ max and as such, letting k = τ −1 r in addition, we obtain We have defined k = τ −1 r , where the constant τ r is the relaxation time, a measure of the relaxation timescale of the chains according the evolution of Γ(τ ). Therefore when τ = τ r , The delay time measures the time by when the polymer has lost its memory of its initial conformation, acting as a timescale by which to measure chain relaxation dynamics. Furthermore, the value of Γ max gives a sense of how free the polymer is to explore different conformations, i.e. it is a measure of deformability. We calculate Γ max for τ = 1.5 ns.

Reparametrization of the F8BT Force Field
Initial structures of 'F0BT' chains (i.e., F8BT chains with their octyl side chains replaced by H atoms) were obtained using the semi-empirical xTB-GFN2 method as implemented in the xTB package. 2 Note that the octyl chains present in F8BT were not included in these calculations since they have been shown to hinder convergence and accuracy of geometry optimisations. 3 Furthermore, they were not included in the F0BT chains used to calculate partial charges for similar reasons. Figure S5: Chemical structure of the F0BT monomer. Note that hydrogen atoms are implicit.

Reparametrization of Partial Charges
F0BT oligomers were used to determine partial charges using the RESP method. Initial structures of F0BT oligomers were minimized using the xTB-GFN2 semi-empirical method before DFT energy minimization was performed. Full geometry optimisations were conducted on each F0BT chain. Using the outputted electron density, the Multiwfn package was used to calculate partial charges using the widely-used restrained electrostatic potential (RESP) method. 4 No constraints were applied when generating the partial charges, with respect to symmetry or with reference to the normal behaviour of the CHARMM force fields. Subsequently, the octyl alkyl chains were parameterised using the CHARMM General Force Field and merged with the calculated charges for the rest of the molecule. Each alkyl carbon atom that is connected to a F0BT carbon had its charge modified simply as q C,new = q C,orignal + q H,removed . Partial charges were then directly implemented in the GRO-MACS topology file, which is available as part of the Supporting Information.
An interesting comparison can be made between the CHARMM36 and RESP partial charges for the F0BT monomer. In donor-acceptor polymers such as F8BT, one can consider the donor unit (here fluorene) to undertake partial charge transfer towards the acceptor unit (here benzothiadiazole). This is reflected in the large, positive partial charge of C9 as calculated by RESP, which is not reflected in the CHARMM partial charge.  Figure S6: F0BT monomer partial charges from RESP and CHARMM. Atom names are the same as in Figure S5. Note that CHARMM underestimates the N partial charges and does not account for partial charge transfer from fluorene to benzothiadiazole, as indicated clearly by the respective C9 partial charges.

Reparametrization of Covalent Bond and Bond Angle Terms
We performed relaxed geometry scans to obtain energy profiles for covalent bond and bonded angle interactions that are poorly described by the classical force field, either to flagged by the large penalties arising in the initial parameterization, or because the chemical analogy made in the initial parameterisation was not appropriate. The geometry scans were performed on an F0BT monomer.
1. Calculate an approximate initial geometry using the semi-empirical GFN2-xTB method with the xTB package 2 2. For each bonded interaction of interest (bond or angle), calculate the potential energy surface using relaxed scans by DFT calculations (B3LYP exchange-correlation functional 5,6 and def2-TZVPP Karlsruhe basis set) using the Orca quantum chemistry package 7 3. Fit the target potential energy surfaces with an appropriate classical force field function (we used non-linear least squares as implemented in the scipy Python package) Covalent bonds are described in the CHARMM force field as V (r) = k r (r − r 0 ) 2 , with k r having units of kcal mol −1Å−2 and r 0 having units ofÅ, which we also used for the fitting procedure. However, direct implementation into the GROMACS topology file requires   Bonded angle deformations are described in the CHARMM force field as V (θ) = k θ (θ − θ 0 ) 2 , which we also used for the fitting procedure. However, for direct implementation into the GROMACS topology file requires V (θ) = k θ 2 (θ − θ 0 ) 2 , with k r having units of kJ mol −1 rad −2 and r 0 having units of deg.  Angle Reparametrization of the Donor-Acceptor Dihedral Angle We followed a similar approach to that of Wildman and co-workers to obtain suitable parameters to accurately describe the donor-acceptor dihedral angle.  Figure S9 shows the 'push-pull' intermonomer dihedral pair used in the dihedral implementation in this work. The DFT and MD background energy scans were performed on the left-hand dihedral (note that this choice is arbitrary). Figure S9: The 'push-pull' intermonomer dihedral pair.
Dihedral rotations are described using the following terms in the force field, with force constants and phase differences listed in Table 3. Figure S10 shows the results from each section of the dihedral angle fitting procedure.
Note the large contribution made by the background effects of the classical FF (gray line) and that the difference between the DFT dihedral energy (blue) and the overall classical FF dihedral energy as implemented in this work (green) is within chemical accuracy (i.e., a difference of less than 4 kJ mol −1 ) for all dihedral angles.

Validation of the F8BT Force Field
While direct comparisons to the experimental structure of the amorphous state are not possible given a lack of experimental methods to probe the structure at the atomistic level of detail, experiments do however provide a suitable point of reference by which to test the validity of the new force field. F8BT is insoluble in water but highly soluble in tetrahydrofuran (THF). 8 To test the solubility of F8BT as modelled by our new force field we performed two simulations, both containing 17 F8BT chains in water and THF respectively.
Simulation Details. F8BT was modeled using the reparameterized CHARMM-based force field 9,10 developed as part of this work. Before use in any production simulation, the F8BT chain was subject to high temperature dynamics in vacuo to randomize its intermonomer dihedral distribution. We note that frequent transitions between different dihedral states are observed in all simulations in addition. Water was treated with the CHARMMmodified TIP3P model 11 and tetrahydrofuran (THF) was modeled using the CHARMM36 force field. 9,10 The MD simulations were performed using the GROMACS 2019 simulation engine. 12,13 In all simulations, the Lennard-Jones and Coulomb interaction cut-off distances were set to 12Å. The particle-mesh Ewald method was used to calculate long-range electrostatic interactions. Periodic boundary conditions were applied in all dimensions in all simulations. A 1 fs timestep was used for all equilibration and production simulations. The leapfrog integration scheme was used in both simulations. For both water and THF simulations, 17 F8BT chains were randomly placed in a simulation box measuring 110×110×110Å 3 .
Solvent molecules were added to yield a system containing F8BT with a concentration of 10.0 wt. % (see Table 4 for details). Both systems were once more subjected to energy minimization by steepest descent before the temperature was equilibrated to 350 K for 100 ps using the Nosé-Hoover thermostat. Subsequently, the production simulation for the F8BT chain in water was performed using the velocity-rescale thermostat 14 (target temperature of 353 K) and the Parrinello-Rahman barostat 15 (target pressure of 1 atm) for 200 ns to allow for any molecular aggregation to be observed.  Figure S11(c).
Summary of Simulations. We find that our force field correctly describes the solubility of F8BT in water and THF. Figure S11

Theoretical Modelling of the Ring Stacking Network
The following is a summary of key concepts that go into the analysis of inter-chain percolation in a glassy configuration of a system of conjugated polymers. The analysis conceives of the amorphous polymer mixture as a network of randomly interlinked polymers, with the links realized through ring stacking interactions as described in the main paper. We will in the following, therefore, refer to polymers also as nodes or vertices of a network or graph, and to the ring stacking interactions as the edges. The theory described below is based on the cavity-approach to percolation formulated by Karrer et al. ,20 supplemented by techniques to expose the full heterogeneity in the problem and to analyse it in the thermodynamic limit as described in Kühn and Rogers. 21

Cavity Approach to Percolation
As input to our analysis, we use the empirically determined probability mass function (pmf) of the polymer stacking degree determined from the simulation, as displayed in Figure 3(b) in the main text. We will assume that the configuration is a random realization of an amorphous configuration of a network giving rise to the pmf of polymer stacking degrees, and that this stacking degree distribution is characteristic of the true limiting stacking degree pmf one would observe in the thermodynamic limit N → ∞ of infinite system size.

Probabilistic Description
Our analysis is based on a probabilistic analysis of bond percolation in large heterogeneous networks, devised in Karrer et al. 20 They consider a network or graph G = (V, E), with V denoting the set of vertices or nodes of the graph, and E the set of edges or links, the polymer stacking interactions in the case at hand. We use N = |V | to denote the size of the graph and i ∈ {1, . . . , N } to label the nodes.
Bond percolation is a process by which each edge (ij) ∈ E is either kept with probability p or deleted with probability 1 − p. A probabilistic description of percolation is obtained by noting that the configuration of connected clusters will be random and be determined by the random realization of the percolation process. For the purpose of the present study, we will eventually be interested in the results of the probabilistic description in the limit p → 1, where all edges of the system are kept. In in which s ∂i = s j j∈∂i denotes the collection of cluster sizes s j that can be reached through the edges (ij) connected to i, and π For the π (i) j (s j ) a set of self-consistency equations is derived, 20 which express the same line of reasoning, viz.
in which the first contribution corresponds to the configuration that the edge (ij) is absent, while the second contribution describes the configuration in which the edge (ij) is present.
These equations are exact only on trees. However they provide excellent approximations for large finitely connected systems, as the one under study, and are known to become asymptotically exact in the thermodynamic limit. [20][21][22] Equations (6) and (7) are most efficiently analysed using probability generating functions or z-transforms corresponding to π i (s) and π (i) j (s), respectively, defined as 20 and For these, Eqs. (6) and (7) entail so these equations can be solved independently for each z.
Noting that G i (1) is the probability that site i belongs to a cluster of any finite size s, one can conclude that is the probability that site i belongs to the giant cluster, and thus gives the average fraction of sites occupied by the giant cluster.
Given a site i is not on the giant cluster, the expected size of the finite cluster to which i belongs is .
Its evaluation requires the z-derivative H (i) j ′ (z) evaluated at z = 1, which is obtained from (11) as Equations (11) (at z = 1) and (15) can efficiently be solved by forward iteration using random initial conditions, and the solutions can then be used to evaluate (12) and thereby the average fractionḡ of sites contained in the giant cluster (13) as well as the average sizes ⟨s i ⟩ of finite clusters (14) containing individual nodes in the system. Averages are here to be understood as averages over many instances of the percolation problem on the same large complex network, so using the present setup avoids simulating many instances of a percolation experiment for the purpose of averaging. Given a large single instance of a complex network, the solution of (11) will in general be heterogeneous across bonds (ij) entailing that the g i as well as the ⟨s i ⟩ will be heterogeneous across nodes. Note, however, that the analysis described above requires complete knowledge of the original graph G = (V, E). In the absence of such knowledge, one can resort to an analysis of the system in the thermodynamic limit N → ∞ that only requires knowledge of the limiting degree distribution p = (p k ) k≥0 of the nodes in the system.

The Large System Limit
Assuming that, in the limit of infinite system size, there is a limiting joint probability density one can use stochastic self-consistency arguments in Eqs. (11) and (15), as explained by Kühn and Rogers, 21 to obtain the following non-linear integral equation for this joint probability density functioñ Here k c p k with c = ⟨k⟩ is the probability that a randomly chosen neighbour of a node has degree k, assuming as we do that there are no degree correlations in the system. From its solution, which can be very efficiently obtained using a stochastic so-called population dynamics algorithm, 23 one obtains, with reference to Eq. (12), the probability density function P (g) of the g i 21 as and, with reference to Eq. (14), the probability density function Q(s) of the ⟨s i ⟩ as Note that the probability density function Q(s) of the average sizes of clusters containing a randomly selected node is related to the cluster size distribution P (s) via which follows from the fact that each cluster of size s has s nodes which belong to it.
For the purposes of the present investigation we are interested in the p → 1 limit of these identities and results. In this limit, the solution of the system Eqs. Eq. (12), will in this limit themselves be indicator variables, taking the value g i = 1 if i is part of the GCC, and g i = 0, if it is not, with Eq. (12) in this limit expressing the logic that for a node to belong to the GCC, it must be connected to the GCC through at least one of its neighbours. In this sense the system Eqs. (11) of self-consistency equations at z = 1 provide a local algorithm that is able to reveal whether a node belongs to the GCC of a given network or not. One can similarly convince oneself that in the limit p → 1 the solutions (15) and (11) (for z = 1) of equations take values on the non-negative integers N 0 = {0, 1, 2, . . . } entailing that both P (s) and Q(s) will be supported on the positive integers. Results for P (s) obtained from the large system limit of the theory as described above are displayed in Figure 3(c) in the main text. The only input that goes into obtaining them is the empirical degree distribution p = (p k ) of the polymer stacking interaction obtained from the MD simulations of amorphous F8BT which is displayed in